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A..  K*  >  Wfc 


FORMULATION  OF  A  MODAL-SPLIT-EXPLICIT  TIME  INTEGRATION  METHOD 

FOR  USE  IN  THE  UCLA  ATMOSPHERIC  GENERAL  CIRCULATION  MODEL 

1.0  INTRODUCTION 

In  the  explicit  time  integration  of  primitive  equations, 
the  time  step  is  limited  by  the  Courant-Friedrichs-Lewy  ( CFL) 
condition  (Courant  et  al.,  1928),  bas^d  on  the  fastest  phase 
speed  of  waves  allowed  in  the  model,  which  is  the  speed  of  the 
external  gravity  wave;  thus  the  integration  is  rather  slow. 

Since  the  meteorologically  significant  waves  have  speeds  much 
slower  than  that  of  the  external  gravity  wave  and  possess  practi¬ 
cally  all  of  the  energy  of  the  atmospheric^motion,  the  use  of  the 
explicit  scheme  is  rather  uneconomical.  Two  major  approaches 
have  been  taken  to  circumvent  this  difficulty.  First,  since  the 
gravity  waves  are  quasi-linear ,  the  portion  of  the  primitive 
equations  that  govern  the  linear  gravity  waves  can  be  integrated 
implicitly,  and  the  rest  of  the  equations  can  be  integrated 
explicitly  with  a  time  step  dictated  by  the  CFL  condition  for 
slow  moving  waves.  Known  as  the  semi-implicit  method,  it  has 
been  used  both  in  grid  point  models  (Kwizak  and  Robert,  1971)  and 
in  spectral  models  (Robert,  1969)  with  significant  savings  in 
computing  time  and  acceptable  accuracy  (Robert,  et  al.,  1972). 
However,  since  the  speed  of  the  gravity  waves  is  reduced  in  the 
implicit  integration,  its  accuracy  is  questionable  in  the  regions 
of  gravity  waves  excitation  (Messinger  and  Arakawa,  1977)  .  The 
nemi-implicit  method  can  be  modified;  the  terms  governing  the 
linear  gravity  waves  can  be  separated  into  different  vertical 
modes.  Only  modes  with  phase  speeds  greater  than  those  of  the 
meteorological  waves — usually  the  external  gravity  wave  and  the 
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first  two  or  three  internal  gravity  waves — need  to  be  integrated 
implicitly  (Burridge,  1978) . 

The  second  approach  to  expedite  the  integration  involves 
separating  the  terms  related  to  the  gravity  waves  from  the  rest 
of  the  equations.  One  group  of  terms  is  integrated  first  with 
an  optimal  time  step,  either  explicitly  or  implicitly,  and  the 
result  of  this  step  is  used  at  the  beginning  of  the  marching  for 
the  other  group,  employing  a  different  optimal  time  step.  This 
is  called  the  splitting  method  (Marchuk,  1974;  Gadd,  1978). 

When  explicit  schemes  are  used,  the  splitting  method  does  not 
have  the  adverse  effect  of  retarding  the  gravity  waves.  However, 
since  the  different  dynamic  processes  are  calculated  one  at  a 
time,  the  truncation  errors  of  all  steps  are  multiplied  rather 
than  added. 

In  a  logical  progression  Madala  (1980)*  proposed  a  modal- 
split-explicit  scheme  (MSES) ,  which  combines  the  advantages  of 
the  aforementioned  methods,  yet  avoids  their  shortcomings.  In 
this  method  terms  in  the  primitive  equations  governing  linear 
gravity  waves  are  decomposed  into  different  vertical  modes. 

Modes  with  phase  speeds  greater  than  the  meteorological  waves 
are  integrated  explicitly,  with  time  steps  allowed  by  the  res¬ 
pective  CFL  conditions  and  are  recombined  at  periodic  intervals. 
The  other  modes  and  the  rest  of  the  equations  are  integrated 
explicitly  with  a  time  step  determined  by  the  speed  of  the  slow 
moving  waves,  which  is  usually  five  times  as  large  as  the  time 
step  allowed  for  the  fastest  moving  waves.  This  constitutes 
great  savings  in  computing  time.  In  addition  since  all  waves 

*Madala  (1981)  used  the  name  "split-explicit;"  the  first  author  added  the 
word  "modal"  to  distinguish  it  from  the  split-explicit  schemes  used  by 
other  authors. 
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are  integrated  explicitly,  no  adverse  effect  on  gravity  wave 
speed  occurs.  Moreover,  since  all  dynamic  terms  are  integrated 
(or  marched)  simultaneously,  the  multiplicative  truncation  error 
characteristic  of  the  splitting  methods  is  avoided. 

Tftus  far,  MSES  has  been  successfully  applied  to  a  tropical 
cyclone  model  (Madala,  1981)  .  The  speed  of  integration  was  in¬ 
creased  by  a  factor  of  three  over  that  using  the  explicit  method; 
this  factor,  however,  is  an  increasing  function  of  horizontal 
resolution.  The  results  also  showed  good  accuracy.  The  present 
paper  presents  a  further  investigation  of  MSES  in  order  to  ascer¬ 
tain  whether  it  is  also  a  viable  method  in  a  global  model,  in 
which  more  complexities  exist,  such  as  high  terrain  and  the 
diminishing  zonal  grid  size  toward  the  poles.  The  model  chosen 
for  this  effort  is  the  1977  version  of  the  UCLA  model  (Arakawa 
and  Lamb,  1977;  hereafter  referred  to  as  AL;  Arakawa  and  Mintz, 
1974)  . 

This  investigation  shows  that  modifications  are  necessary 
(Section  2.2);  that  additional  features  should  be  added  (Section 
2.3);  and  that  the  pressure  averaging  technique  can  be  used  to 
improve  the  efficiency  of  MSES  (Section  2.5).  In  addition,  the 
results  of  a  72-hour  test  run  showed  that  MSES  is  viable,  effi¬ 
cient,  and  accurate  for  global  weather  prediction  purposes. 

2 . 0  FORMULATION 

2 . 1  Basic  Equations 

Since  MSES  treats  the  terms  governing  linear  gravity 


waves  in  the  equations  of  motion  differently  from  the  rest, 
these  terms  must  first  be  identified.  The  vertical  coordinates 


in  the  UCLA  GCM  is 


TT 


where  P  is  pressure;  PT  is  a  constant  pressure  at  the  model  top 
(currently  PT  =  50  mb)  ;  tt  «  Pg  -  PT;  and  Pg  is  the  pressure  at 
the  surface.  Because  there  is  only  one  layer  in  the  strato¬ 
sphere,  the  original  definition  of  o  in  the  stratosphere  (AL, 
p.  207)  is  not  used  for  reasons  of  simplicity.  All  notations 
follow  those  of  AL  unless  otherwise  indicated. 

The  zonal  momentum  equation  in  the  flux  form  ia 

~  Ciru)  +  ^  (ttu2)  +  — ~~  (ituvcos2!)) 

cos  $ 

+  (iui)  -  »fv  -  ■"(!$) 

-  Circa)  +  HFX  ,  (2) 

where  u  is  zonal  velocity,  v,  meridional  velocity;  <J>,  geopoten¬ 
tial;  a,  specific  volume;  f,  Coriolis  parameter;  <J>,  latitude, 

X,  longitude;  a,  earth's  radius;  F  ,  zonal  frictional  force; 
o  =  da/dt;  9x  =  acos<}>3X;  and  3y  =  a3<j>. 

Since  0  =  4>  +  <J>  ,  where  cj>  is  the  surface  geopotential 

S  a  S 

and  <J>  is  the  <t>  measured  from  surface,  the  pressure  gradient 

a 

force  terms  in  Eq.  (2)  can  be  written  as 

-  TZ  t,rV  +  [cv5™*  +  lVow)']  A  ”  -  *  TZ 


where  the  double  bar  denotes  global  average  on  a  o-surface,  and 


the  prime  denotes  the  deviation  from  it.  If  <f>  =  n <j>  - 

- —  a 

(<j>  -cn ra)Tr,  Eq.  (2)  becomes 

cl 

It  (7TU)  +  I?  I  =  V  (3) 

where  Ru  denotes  remaining  terms,  which  vary  slowly  in  time  com¬ 
pared  to  the  terms  on  the  left  side  of  Equation  (3).  The  addi¬ 
tional  operator,  denoted  by  the  dash  bar,  suppresses  the  ampli¬ 
tude  of  the  short  waves  to  overcome  the  problem  associated  with 
the  diminishing  zonal  grid  size.  This  operator  is  fully  des¬ 
cribed  in  AL,  p.  248,  and  should  be  applied  to  the  zonal  pressure 
gradient  term  and  to  the  zonal  mass  flux  term.  In  a  similar  man¬ 
ner,  the  equation  governing  meridional  momentum  can  be  written 
as : 


9  ,  ,  ,  9  —  D 

9t  (lTv)  +  3y  ♦  =  V 

The  surface  pressure  tendency  equation 
|^T  7T  +  ^N2>T  <D>  =  0, 


is. 


(4) 


(5) 


where  D  =  V*  (nv)  =  (  ttu  )  +  (ttvcos4>)  (AL,  Eq.  (166)). 

T 

The  row  vector  <N2->  is  given  in  the  appendix. 

The  thermodynamic  equation  (AL,  Eq.  (209))  is 

It  (7rV  +  V*<7rVkV  +  [(7rd)k+b  Pk®k+is 

(7r<J)k-Js  ^k-,]  =  (ona)k  +  vk-v)  tt  +  J-  Qk , 

(6) 
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where  the  subscripts  are  level  indices;  6,  potential  temperature; 
P  K 

p  s  T/e  s  (lO'ffff)  '  and  cp'  sPeci^ic  heat  at  constant  pressure. 
Integration  of  the  continuity  equation  gives  <tto>  =  [N^]  <D> . 
Thus,  Eq.  (6)  can  be  written  as: 


or  as: 


^  U<T>)  +  tM2]  <D>  »  <Rt> 


(7) 


where  [M^]  <D>  denotes  the  linear  part  of  the  last  two  terms  on 
the  left  side  of  Equation  (6) ,  and  RT  denotes  all  the  remaining 
terms.  The  matrices  [M2]  and  [M4]  are  given  in  the  appendix. 
Instead  of  n ,  tt  is  used  in  computing  [M4]. 

Equations  (3),  (4),  (5),  (6)  and  the  hydrostatic  equation 
form  the  complete  set  of  governing  equations  in  the  model.  When 
R  =  Rv  =  Itj,  =  0 ,  these  equations  govern  the  linear  gravity 
waves  in  an  atmosphere  with  no  basic  motion. 


2.2  The  Modal-Split-Explicit  Method 

In  the  UCLA  GCM,  $  is  related  to  the  temperature  nonlin- 
early  (AL,  Eqs.  (207) -(208) ) ,  and  this  will  be  denoted  by  a  sub¬ 
script  N,  thus  (<fe  ).T.  However,  one  important  requirement  in  the 
present  scheme  is  that  <j>  must  be  linearly  related  to  tempera- 

a 

ture.  This  requirement  is  met  by  a  previous  definition  of  <J>  in 

a 

the  UCLA  model  (Arakawa,  1972;  and  Appendix) ,  and  this  will  be 
denoted  by  a  subscript  L,  thus  (4>a)L«  Accordingly,  the 
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hydrostatic  equation  has  the  form:  <if  >T  =  [M.]  <T> ,  where 

a  l,  -L 

[M  ]  is  given  in  the  appendix.  The  matrix  is  a  function  of 

tt  ,  when  PT  in  Eq.  (1)  is  not  equal  to  zero;  therefore,  it  is  a 
function  of  the  horizontal  coordinates.  As  part  of  the  simpli¬ 
fication,  [M^]  computed  with  u  will  be  used  for  all  locations. 

In  this  way,  the  definition  of  <<p  >  is  modified. 

a  i-j 

Following  the  incorporation  of  the  foregoing  arguments 
Equation  (3)  should  be  rewritten  as: 


5 

3 1 


(r.u) 


R  . 


u 


(8) 


When  using  leapfrog  scheme  in  time,  the  finite  difference  form 
of  (3)  over  a  2Ax  interval  can  be  written  as, 


(TTu)t-At+{n  +  1)AT 


(nu) 


t-At+ (n-1) Ax 


+  2  Ax 


3 _ 

9x 


---t-At+nAx 


=  2  Ax  K 


t-At+nAx 


(9) 


where  Ax  is  a  time  interval  allowed  by  the  CFL  condition  based 
on  the  speed  of  fastest  waves,  and  At  is  that  based  on  the 
meteorological  waves.  The  subscript  N  indicates  that  it  is  com¬ 
puted  using  nonlinear  form  of  hydrostatic  equation.  Marching 
Eq.  (9)  with  2 Ax  intervals  between  t-At  and  t+At  and  averaging 
these  equations  in  time  given, 


Ittu) 


t+At 


(7Tu)t  At  +  2At 


3 

3x 


— 2-2At 

(  ±  >N 


2At  R^At 


=  2 At  R 


u 


(10) 
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9  A  r 

where  (~)  is  the  discrete  averaging  over  the  2At  interval. 

With  regard  to  the  stability  of  the  scheme,  the  approxima- 
%  2At  t 

tion  of  Ru  by  Ru  is  allowed,  since  Ru  does  not  create  a 
high  frequency  variation  in  ttu,  and  since  it  is  a  slow  vary¬ 
ing  term.  This  approximation  introduces  a  truncation  error, 

which  is,  however,  bearable.  When  the  term  2At  JLf*2At_52At  -rt\ 

3x  L  -N  — L/ 

is  added  to  both  sides,  Eq.  (10)  becomes 

Uu)t+At  -  (nu)t-At  +  2At  ±  (|A€-g) 

2At  Ru  +  2At  ^  -iN  -ij 

24t  Eu  +  2At  4 


2At  Ru  -  24t  H  4 


(11) 


The  right  side  of  Eq.  (11)  gives  the  change  of  ttu  over  the 
2At  interval  in  one  leapfrog  step.  Therefore,  Eq.  (11) 
becomes : 

(7TU)t+At  -  (7TU)t_At  +  2A  t  (gAt4^ 


=  (ttu) 


t+At 

E 


(ttu) 


t-At 


(12) 


where  the  subscript  E  denotes  the  results  from  marching  once  with 

the  2At  interval  without  the  help  of  the  MSES  scheme.  Thus, 

(nu) t+At  can  be  found  by  marching  Eq.  (8)  over  2At  only  once,  if 
At  . 

<1>L  is  known. 
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Similarly,  the  other  governing  equations  are: 


(rrv)t+At  -  ( nv)  +  2ht 


,  . t+At  ,  .  t-At 

=  Cttv)e  ~  <"v) 


<iTT>t+At  -  <TTT>t_At  +  2 At  [M2]<D2At-Dt> 


=  <irT>^+At  -  <rrT>t  At 
E 


TTt+At  -  TTt_At  +  2 At  <N2>T  <D2At-Dfc> 


t+At  t-At 
=  tte 


Since  the  aim  now  is  to  obtain  the  two  sets  of 
—  2At  ^  2At 

unknowns,  <£>  and  <D>  ,  the  four  sets  of  equations  can 

be  reduced  to  two  sets.  Forming  [M,  J  (14)  +  <ava-<}>*  >  (15) 

j.  L 

gives : 


<i^t+At  “  <i^,t  At  +  2At  [-M-3]<D2At-Dt> 


-X  st+At  -t-At 

<ir.>E  "  <4i.> 
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where  [H^J  =  IMj]  +  <otto.— (}>L>  <Nj>  .  The  grid  scheme 

C  (AL,  p.  182)  used  in  the  UCLA  model  is  not:  changed. 

In  order  to  form  the  divergence  equation  from  Egs. 
(12)  and  (13)  it  is  convenient  to  define  ir  at  the  u  (v)  points 
as  the  average  of  the  two  neighboring  it’s  in  the  zonal 
(meridional)  direction.  However,  it  should  be  emphasized 


that  the  original  definition  of  ir  at  the  u  and  v  points  (AL, 
p.  242)  are  still  used  in  computing  and  v^+At. 

Equations  (12)  and  (13)  give  the  divergence 


equation 


<D>t+At  -  <D>t“At  +  2Lt  V2  <<J>2it-<J>^> 

L  — 1* 


-  <D>t_At  , 


(17) 


where  V2  =  (cos <J>  Note  that  the  dash 

bar  operator  is  used  twice.  Spectral  equations  governing  the 
amplitudes  of  the  natural  gravity  modes  of  the  numerical  model 
can  be  obtained  by  premultiplying  equations  (16)  ard  (17)  by 
[E]~  ,  where  [E]  is  a  matrix  whose  columns  are  the  eigenvectors 


and 


(<DE>^T_<DE>t)  _  (<DE>T-AT_(DE)tj 


O  _ p  1  _p  u 

+  2Atv  -4)^  > 


=  (AT/At)[E]-1(<D>‘+it-<D>t-it)  ,  (21) 


where  the  right  sides  of  these  equations  are  not  changed 
(except  for  the  Ax/At  factor),  since  they  vary  rather  slowly 
over  a  2At  period. 

Equations  (20)  and  (21)  are  marched  in  the  2At 

interval  with  different  Ax  for  each  mode  as  determined  by 

the  CFL  condition  based  on  the  phase  speed  of  that  mode. 

The  right  sides  of  Eqs.  (19)  and  (20)  are,  of  course,  held 

=2At  Et 

constant  in  this  marching.  The  quantities  (D^  -D  )  and 
-E2At  -Efc 

(£  )  are  obtained  from  averaging  the  results  of  this 

marching.  Only  modes  with  phase  speeds  greater  than  the 
maximum  speed  of  the  meteorological  waves  have  to  be  marched 

i2At  Et 

in  this  manner.  Finally  E  is  multiplied  with  (D  -D  ) 

and  (j[E  -£E  )  to  obtain  (D2At-Dt)  and  (£2  t-^t)  respectively. 
L  L  L  L 

In  summary,  the  procedure  of  the  MSES  method  is: 

(1)  Calculate  the  matrices,  and  the  eigenvalues 
and  eigenvectors  of  IM3I,  using  the  prede¬ 
termined  global  mean  quantities.  This  step 
is  done  only  once. 
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I  _ p  t~At  _ _  t\  /  E  p 

(2}  Compute  “<iL>  /and(<D  *  -<D  >  J. 

C3)  Use  the  larger  time  interval  At,  march 

forward  one  step  and  compute  the  right  hand 
sides  of  Eqs.  (20)  and  (21). 

(4)  March  Eqs.  (20)  and  (21)  with  different  At 
for  different  modes.  For  the  six-level 

UCLA  model  (AL,  p.  176) ,  only  the  first  three 
modes  need  to  be  integrated  with  Ax  =  (1/5, 
1/3,  1/2) At. 

(5)  Average  the  result  from  step  (4)  to  obtain 
(<I>LAt'<i>i,)  and  (<3>2At“<D>t)  which  then 

allow  the  calculation  of  the  predicted  quan¬ 
tities  at  t+At  in  Eqs.  (12)  throuqh  (15). 

Although  the  above  descriptions  are  based  on  the 
leapfrog  scheme,  the  changes  are  minor  when  using  the  Matsuno 
scheme,  which  is  periodically  used  in  the  UCLA  GCM,  (AL,  p. 
260).  In  the  first  part  of  the  Matsuno  scheme,  the  super¬ 
scripts  t+At,  t-At  and  2At  in  the  above  equations  are  changed 
to  t*,  t  and  At,  respectively.  Then,  in  the  second  part  of 
the  Matsuno  scheme,  the  superscripts  t-At  and  2At  are  changed 
to  t*  and  At,  respectively.  In  both  parts,  of  course,  the 
factor  2At  is  changed  to  At. 

2. 3  Modification  of  the  Advective  Terms  Near  the  Poles 

Even  when  At  is  increased  by  a  factor  of  five,  the 
dash  bar  operator,  described  in  subsection  2.2,  is  sufficient 
to  circumvent  the  problem  of  linear  instability  related  to 
the  gravity  wave  terms,  due  to  the  diminishing  zonal  grid 
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toward  the  poles.  However,  the  advective  terms  3/3x  jOFi7)gj 
— where  q  denotes  u,  v  and  T--can  create  linear  instability 
near  the  poles  when  At  is  enlarged.  The  time  increment 
allowed  for  the  advective  terms.  At,  is  limited  by  At  <  eAx/ 
(u+c^) ,  where  cm  is  the  speed  of  the  meteorological  waves, 
u  is  that  of  the  basic  flow,  and  e  has  a  magnitude  of  the 
order  of  1,  with  its  exact  value  depending  on  the  time  dif¬ 
ferencing  scheme  used,  when  u  is  large  in  the  polar  region, 
where  Ax  is  small.  At  cannot  be  increased  by  a  factor  of  five. 
One  obvious,  yet  undesirable,  solution  is  to  use  a  smaller 
At.  The  solution  used  here  is  to  apply  the  dash  bar  opera¬ 
tor  to  q  in  the  terms  8/3x  [(ifu)q]  in  Eqs.  (272)  and  (299) 
of  AL.  Thus,  in  step  2  of  the  MSES  procedure  in  subsection 
2.2,  the  second  term  in  Eq.  (299)  is  changed  to  (with  notation 
of  AL) : 

- ? 

[F (  T  )  ]. 

Also,  the  changes  in  Eq.  (272)  are  to  substitute  Eq.  (275) 
into  Eq.  (272)  and  then  to  apply  the  dash  bar  operator  on 
u  in  those  terms  containing  F*.  Similar  changes  are  made 
for  the  v  component  equation.  Short-term  tests  (72  hrs.) 
show  this  method  is  successful  in  controlling  linear  insta¬ 
bility;  however,  since  the  quadratic  conservation  properties 
no  longer  hold,  there  can  be  some  side  effects  in  this 
approach,  as  will  be  shown  in  the  result  section. 
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2. 


The  Time  Differencing  Method 


As  shown  in  Figure  1,  the  time  differencing  method 
used  to  calculate  the  dynamic  terms  consists  of  a  series  of 
the  leapfrog  schemes,  with  a  periodic  insertion  of  the 
Matsuno  scheme.  The  source  and  sink  terms,  and  the  verti¬ 
cal  flux  convergence  term  of  the  moisture  equation,  are 
calculated  as  an  instantaneous  adjustment.  For  At  =  30  min., 
these  calculations  are  done  after  each  leapfrog  or  Matsuno 
step.  Thus,  the  time  differencing  method  for  the  physical 
processes  is  not  changed. 


(4,T+1_<t>T_1)/2At  +  gH  ux  =  0  , 

and 

(uT"1-uT"1)/2At  +  ^  =  0  , 

where 

(V  =  o[(  )T  +  1  +  {  )T~1\  +  (l-2a)  (  )T  . 

with  this  method,  At  can  be  almost  twice  as  large  as  that 
allowed  by  the  CFL  condition.  Thus,  when  this  method  is 

p 

applied  to  Eqs.  (20)  and  (21)  the  term  £  is  changed  to 
£  .  In  the  present  investigation,  a  =  0.24  is  used.  Thus, 

Eqs.  (20)  and  (21)  can  be  integrated  with  Ax  =  (1/3 , 1/2 ,1/1 . 5) At 
for  the  first  three  modes.  A  detailed  analysis  of  this 
method  was  performed  by  Brown  and  Campana  (1978). 


3.0  RESULTS  AND  DISCUSSION 

The  performance  of  the  present  MSES  formulation 
is  of  course  judged  by  its  stability,  efficiency,  and 
accuracy.  The  following  discussions  are  directed  toward 
these  characteristics. 

3 . 1  Stability 

Two  12-hr.  and  one  72-hr.  runs  starting  from 
different  initial  conditions  showed  that  the  present  MSES 
formulation  is  capable  of  avoiding  linear  instability. 
Since  the  linear  stability  criterion  depends  not  only  on 
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the  wave  speeds  but  also  on  the  wind  speed,  present  formulation 
should  be  tested  with  large  wind  speed,  especially  in  the  polar 
region  where  zonal  grid  size  is  small.  In  one  preliminary  12-hr. 
test  run  the  wind  speed  near  the  nortnpole  reached  more  than 
50m/s  and  no  linear  instability  occurred. 

3 . 2  Efficiency 

Tests  show  that  when  the  UCLA  GCM  with  A <p  =  4°, 

AA  =  5°  and  At  =  30  min.  (five  times  of  the  original  At)  was  run 
on  a  CDC  175  without  physical  processes,  the  reduction  in  inte¬ 
gration  CPU  time  by  MSES  is  48%  of  which  4%  is  due  to  the 
addition  of  the  pressure  averaging  methods.  The  efficiency  can 
be  greater  if  the  resolution  increases.  The  net  reduction  in 
CPU  time  for  the  full  model  depends  on  the  amount  of  the  CPU  time 
used  for  the  physical  processes,  which  of  course  get  no  improve¬ 
ment  from  MSES.  The  reduction  is  24%  for  the  full  UCLA  GCM  with 
the  aforementioned  grid  size.  Only  when  a  large  portion  of  CPU 
time  is  allocated  to  dynamic  processes  does  it  pay  to  adopt  the 
MSES. 

3 . 3  Accuracy 

Two  72-hr.  runs  with  and  without  MSES  were  made 
starting  from  the  same  well-balanced  initial  conditions,  which 
is  the  end  of  another  72-hr.  run.  The  initial  conditions  and 
the  results  at  hour  72  are  shown  in  Figures  2  and  3  for  sea 
level  pressure  and  50U  mb  geopotential  height,  respectively .  The 
difference  between  the  two  runs  is  very  small  when  compared  with 
the  changes  over  72  hrs.  The  RMS  differences  in  surface  pres¬ 
sure,  in  500  mb  geopotential  height,  and  in  500  mb  zonal  wind 
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between  the  two  runs  are  given  in  Table  1.  The  RMS  difference 
in  500  mb  geopotential  height  at  hour  72  is  16  meters,  which  is 
very  small,  when  compared  with  a  typical  forecast  error  of  75 
meters.  Also  the  RMS  difference  in  500  mb  zonal  wind  at  hour  72  is 
2.4  m/s  compared  with  the  typical  48  hr.  forecast  error  of  8.5 
m/s.  Although  there  was  concern  whether  computing  the  matrices 
using  a  globally  averaged  n  was  appropriate  for  high  terrain 
areas,  no  adverse  effect  over  these  regions  was  detected.  One 
obvious  discrepancy  between  the  two  runs  is  that  after  72  hrs., 
a  strong  high  centered  at  35°  E,  80°  S,  appeared  in  the  run 
using  MSES,  whereas  only  a  weak  high  occupied  this  region  in  the 
other  run.  A  separate  72-hr.  run  using  MSES,  and  starting  from 
the  same  initial  conditions  but  without  the  physical  processes, 
also  has  this  intense  high  (Figure  4) .  Thus,  it  is  reasonable 
to  consider  this  discrepancy  a  side  effect  of  the  modification 
described  in  Section  2.3. 

Table  1.  RMS  Differences 


RMS  Difference 

Surface 
Pressure  (mb) 

500  mb 

Geopotential 
Height  (m) 

500  mb 
Zonal 

Wind  (m/s) 

24 

hr . 

1.09 

8.93 

1.03 

48 

hr. 

1.54 

12.40 

1.79 

72 

hr . 

1.97 

15.78 

2.40 
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Overall,  the  results  are  remarkably  good  outside  the 
regions  very  close  to  the  poles.  The  poor  performance  near  the 
poles,  however,  is  not  a  direct  consequence  of  applying  MSES, 
but  it  is  related  to  the  horizontal  differencing  scheme,  and  to 
the  diminishing  zonal  grid  size  near  the  poles.  Whether  this 
problem  precludes  long-term  integration  was  not  answered  in  this 
study,  due  to  the  limitations  on  computing  resources. 

Thus  far  the  discussion  on  accuracy  is  based  on  the 
runs  starting  from  a  well-balanced  initial  condition.  Whether 
the  good  accuracy  can  still  be  maintained  in  runs  starting  from 
an  initial  condition  that  is  generated  from  an  initialization 
program  is  again  not  answered  due  to  the  limitations  on  computing 
resources.  Nevertheless,  since  a  reasonably  good  initialization 
program  should  generate  a  fairly  balanced  initial  condition,  this 
problem  may  not  be  a  justified  concern.  (However,  incompatabi- 
lity  between  the  initial  thermodynamic  fields  and  the  cumulus 
parameterization  scheme  can  generate  initial  imbalance.) 

Besides,  should  this  problem  be  serious,  the  use  of  MSES  can  be 
delayed  fora  while  (say,  12  hrs.),  after  the  initial  imbalance 


has  died  down. 
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Fig.  2d  —  72  hr.  sea  level  pressure  (mh),  Southern  Hemisphere,  no  MSES 


Fig.  3a  —  Initial  500  mb  geopotential  (60  meter  interval),  Northern  Hemisphere 
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Fig.  3c  —  72  hr.  500  mb  geopotential  (60  meter  interval).  Northern  Hemisphere,  no  MSES 
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Fig.  3e  —  72  hr.  500  mb  geopotential  (60  meter  interval),  Northern  Hemisphere,  with  MSES 


1 
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Fig.  3f  —  72  hr.  500  mb  geopotential  (60  meter  interval),  Southern  Hemisphere,  with  MSES 


Fig.  4a  —  72  hr.  sea  level  pressure  (mb),  Northern  Hemisphere,  with  MSES,  dry  model 
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APPENDIX 


In  a  K-level  o  coordinate  model,  where 


a  =  CP-PT)/CPS-PT)  , 


<N2>  =  (Ac1,Ao2, - Aok)  , 


The  hydrostatic  equation  for  6  (Arakawa,  1972, 

a 


p.  25)  is: 


UJ  ~  (0a) 

a  k  a  k+1 


CP  l(= 

T 


|(^+i)K  -  (pk)K] 


Tk  +  Vi 


k+1 


where  k  =  R/cp,  and 


(A-l) 


K 


K 


E  ^k  f  Aok  ■  (ok+iBk+ok-iak)  Tk,(A-2) 


k=l 


>  [(¥)'  •  ■ 


for  k  <  K  -  1 
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for  k  =  K  , 


V 


and 


0 


f  or  k 


1 


1 

2 


-(M 


for  k  >  2 


The  matrix  can  be  constructed  from  the  preceding  equa¬ 

tions,  where  Equations  (A-l)  and  (A-2)  give  <f>  R  ( k  -  1,....K) 

as  a  linear  function  of  T,  (k  =  1, - K)  .  The  element  (M^ifj 

equals  the  4>  i ,  computed  using  Tj  =  1  and  Tk  =  0  for  all  k  /  j  • 


[Np  - 


°2^°l  °^°2 

E  E 

a  2  Ao  ^  cr 2 


E, 

-o2^oK 

Ea 


VK+1AC1  0K+lA°2"""°K+iaaKi 


(A-3) 


aE  denotes  a  at  an  interface  between  two  layers  as  shown 
in  Fig.  A-l.  When  <tto>  =  [N1J<D>  is  substituted  in  the 
thermodynamic  equation  M2  and  can  easily  be  construc¬ 

ted  . 

The  thermodynamic  equation  is 


It  <"  V  +  "•>"  Vk>  +  Ao7  "‘“'k+^kV* 


(A-4) 
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1 


where  V 


=  [NjJ 


and 


[m2] 


[M4]  +  [T]  +  i 

P 


It  should  be  emphasized  here  that  all  quantities  in  the 
matrices  are  computed  using  the  globally  averaged  temperature 
and  surface  pressure. 
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